Cel

Celem projektu jest określenie jakie mogą być główne przyczyny stopniowego zmniejszania się długości śledzi oceanicznych wyławianych w Europie.

Zbiór danych zostanie wczytany z pliku CSV, następnie musi zostać poddany wstępnemu oczyszczaniu. # Opis zbioru danych Analiza dotyczy zbióru danych na temat połowu śledzia oceanicznego w Europie. Do analizy zebrano pomiary śledzi i warunków w jakich żyją z ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.

Poniżej znajdują się szczegółowe opisy konkretnych atrybutów:

Nazwa kolumny Opis Dodatkowa Informacja
length długość złowionego śledzia [cm]
cfin1 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1]
cfin2 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2]
chel1 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1]
chel2 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]
lcop1 dostępność planktonu [zagęszczenie widłonogów gat. 1]
lcop2 dostępność planktonu [zagęszczenie widłonogów gat. 2]
fbar natężenie połowów w regionie [ułamek pozostawionego narybku]
recr roczny narybek [liczba śledzi]
cumf łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku]
totaln łączna liczba ryb złowionych w ramach połowu [liczba śledzi]
sst temperatura przy powierzchni wody [°C]
sal poziom zasolenia wody [Knudsen ppt]
xmonth miesiąc połowu [numer miesiąca]
nao oscylacja północnoatlantycka [mb]

Wykorzystane biblioteki

library(knitr)
library(ggplot2)
library(polycor)
library(heatmaply)
library(tidyr)
library(plotly)
library(VIM)
library(caret)
library(klaR)
library(dplyr)

Powtwarzalne wyniki

set.seed(23)

Oczyszczenie danych

Wczytanie danych z pliku CSV

raw_data <- read.csv(file= "sledzie.csv", header= TRUE, sep= ",", na.strings= "?")

Podstawowa analiza arybutów

str(raw_data)
## 'data.frame':    52582 obs. of  16 variables:
##  $ X     : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ length: num  23 22.5 25 25.5 24 22 24 23.5 22.5 22.5 ...
##  $ cfin1 : num  0.0278 0.0278 0.0278 0.0278 0.0278 ...
##  $ cfin2 : num  0.278 0.278 0.278 0.278 0.278 ...
##  $ chel1 : num  2.47 2.47 2.47 2.47 2.47 ...
##  $ chel2 : num  NA 21.4 21.4 21.4 21.4 ...
##  $ lcop1 : num  2.55 2.55 2.55 2.55 2.55 ...
##  $ lcop2 : num  26.4 26.4 26.4 26.4 26.4 ...
##  $ fbar  : num  0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
##  $ recr  : int  482831 482831 482831 482831 482831 482831 482831 482831 482831 482831 ...
##  $ cumf  : num  0.306 0.306 0.306 0.306 0.306 ...
##  $ totaln: num  267381 267381 267381 267381 267381 ...
##  $ sst   : num  14.3 14.3 14.3 14.3 14.3 ...
##  $ sal   : num  35.5 35.5 35.5 35.5 35.5 ...
##  $ xmonth: int  7 7 7 7 7 7 7 7 7 7 ...
##  $ nao   : num  2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...

Zbiór zawiera 52582 rekordów rozmieszczonych w 16 kolumnach (z czego jedna jest kolumną porządkową).

Zmiana kolumn

Zmienna xmonth, która reprezentuje miesiąc połowu powinna zostać zamieniona z ciągłej na kategoryczną, by nie traktować jej jako liczbę. Zmienna totaln, która reprezentuje łączną liczbę ryb złowionych, powinna zostać zmieniona na całkowitą.

raw_data <- raw_data %>% 
  mutate(xmonth= as.factor(xmonth), totaln= round(totaln), totaln= as.integer(totaln))

Problem pustych danych

#Puste dane liczbowo na kolumnę
apply(raw_data, 2, function(x){ sum(is.na(x)) })
##      X length  cfin1  cfin2  chel1  chel2  lcop1  lcop2   fbar   recr 
##      0      0   1581   1536   1555   1556   1653   1591      0      0 
##   cumf totaln    sst    sal xmonth    nao 
##      0      0   1584      0      0      0
#Puste dane procentowo na kolumnę
apply(raw_data, 2, function(x){ sum(is.na(x)) / length(x) })
##          X     length      cfin1      cfin2      chel1      chel2 
## 0.00000000 0.00000000 0.03006732 0.02921152 0.02957286 0.02959188 
##      lcop1      lcop2       fbar       recr       cumf     totaln 
## 0.03143661 0.03025750 0.00000000 0.00000000 0.00000000 0.00000000 
##        sst        sal     xmonth        nao 
## 0.03012438 0.00000000 0.00000000 0.00000000

Zbiór zawiera również wartości puste - te pojawiają się głównie w kolumnach z informacją o dostępności planktonu oraz temperaturze przy powierzchni wody.

aggr(raw_data, plot= TRUE, 
     col= c('#fa9fb5', '#2b8cbe'), 
     numbers= TRUE, 
     prop= FALSE, 
     bars= FALSE, 
     labels= names(raw_data), 
     cex.axis= 0.8, 
     ylab=c("Histogram brakujących danych","Wzorzec"))

Jak zobrazowano na wykresie powyżej, rozkład wartości pustych w kolumnach:

  • cfin1 - dostępność planktonu - skupisko Calanus finmarchicus gat. 1
  • cfin2 - dostępność planktonu - skupisko Calanus finmarchicus gat. 2
  • chel1 - dostępność planktonu - skupisko Calanus helgolandicus gat. 1
  • chel2 - dostępność planktonu - skupisko Calanus helgolandicus gat. 2
  • lcop1 - dostępność planktonu - skupisko widłonogów gat. 1
  • lcop2 - dostępność planktonu - skupisko widłonogów gat. 2
  • sst - temperatura przy powierzchni wody stopnie °C

Usuwając wiersze z wartością NA, utracilibyśby stosunkowo dużo danych - lepszym pomysłem jest zastąpienie wartości brakującej średnią z konkretnego połowu. Bazując na fakcie, iż kolumny totaln, xmonth, nao definiują konkretny połów oraz nie zawierają one żadnych wartości pustych, posłużą one do grupowania. Dane zostały zgrupowane względem połowów, a następnie wartości puste zostały zamienione na średnią z tych połowów.

data <- raw_data %>%
  group_by(totaln, xmonth, nao) %>%
  mutate_each(funs(replace(., which(is.na(.)),
                                mean(., na.rm=TRUE))))

Problem duplikatów

no_x <- data %>% select(-X)
sum(duplicated(no_x))
## [1] 45694

Możemy zaobserwować, iż 45694 rekordów to duplikaty. Pojawiają się one wewnątrz jednego połowu, dlatego usunięcie ich nie wpłynie negatywnie, ani nie sfałszuje danych. Dla uproszczenia grafów oraz dalszych obliczeń, wszystkie duplikaty zostały usunięte, tym samym zbiór danych uszczuplił się do 6888 rekordów.

w_duplicates <- unique(no_x[, 1:15])
w_duplicates <- w_duplicates %>% 
  ungroup %>%
  mutate(X = seq_len(n())) %>% 
  select(X, everything())

Statystyki zbioru danych

Podstawowa analiza arybutów

str(w_duplicates)
## Classes 'tbl_df', 'tbl' and 'data.frame':    6888 obs. of  16 variables:
##  $ X     : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ length: num  23 22.5 25 25.5 24 22 23.5 22.5 22 24.5 ...
##  $ cfin1 : num  0.0278 0.0278 0.0278 0.0278 0.0278 ...
##  $ cfin2 : num  0.278 0.278 0.278 0.278 0.278 ...
##  $ chel1 : num  2.47 2.47 2.47 2.47 2.47 ...
##  $ chel2 : num  21.4 21.4 21.4 21.4 21.4 ...
##  $ lcop1 : num  2.55 2.55 2.55 2.55 2.55 ...
##  $ lcop2 : num  26.4 26.4 26.4 26.4 26.4 ...
##  $ fbar  : num  0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
##  $ recr  : num  482831 482831 482831 482831 482831 ...
##  $ cumf  : num  0.306 0.306 0.306 0.306 0.306 ...
##  $ totaln: int  267381 267381 267381 267381 267381 267381 267381 267381 267381 267381 ...
##  $ sst   : num  14.3 14.3 14.3 14.3 14.3 ...
##  $ sal   : num  35.5 35.5 35.5 35.5 35.5 ...
##  $ xmonth: Factor w/ 12 levels "1","2","3","4",..: 7 7 7 7 7 7 7 6 6 6 ...
##  $ nao   : num  2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...

Zbiór danych po oczyszczaniu zmniejszył się do 6888 wierszy, liczba kolumn pozostała niezmieniona i wynosi 16. ## Analiza arybutów

knitr::kable(summary(w_duplicates))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
Min. : 1 Min. :19.00 Min. : 0.00000 Min. : 0.0000 Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849 Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137 Min. :12.77 Min. :35.40 8 : 879 Min. :-4.89000
1st Qu.:1723 1st Qu.:24.00 1st Qu.: 0.02778 1st Qu.: 0.2500 1st Qu.: 2.469 1st Qu.:15.030 1st Qu.: 2.5479 1st Qu.:20.094 1st Qu.:0.1580 1st Qu.: 364794 1st Qu.:0.11008 1st Qu.: 307276 1st Qu.:13.64 1st Qu.:35.51 10 : 879 1st Qu.:-1.69000
Median :3444 Median :25.50 Median : 0.14158 Median : 0.3714 Median : 4.811 Median :21.435 Median : 5.9167 Median :24.859 Median :0.3320 Median : 459347 Median :0.21476 Median : 539558 Median :13.98 Median :35.51 7 : 747 Median : 0.20000
Mean :3444 Mean :25.32 Mean : 0.55913 Mean : 1.7403 Mean : 8.801 Mean :21.157 Mean : 11.3557 Mean :27.683 Mean :0.3202 Mean : 543028 Mean :0.21417 Mean : 523418 Mean :13.94 Mean :35.52 9 : 680 Mean : 0.08938
3rd Qu.:5166 3rd Qu.:26.50 3rd Qu.: 0.36032 3rd Qu.: 1.5701 3rd Qu.: 9.667 3rd Qu.:26.324 3rd Qu.: 12.4959 3rd Qu.:35.153 3rd Qu.:0.4250 3rd Qu.: 774993 3rd Qu.:0.28116 3rd Qu.: 763083 3rd Qu.:14.21 3rd Qu.:35.52 6 : 559 3rd Qu.: 1.80000
Max. :6888 Max. :32.50 Max. :37.66667 Max. :19.3958 Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736 Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595 Max. :14.73 Max. :35.61 5 : 521 Max. : 5.08000
NA NA NA NA NA NA NA NA NA NA NA NA NA NA (Other):2623 NA

Rozklad wartosci

data_dist <- w_duplicates %>% 
  select(-X) %>% 
  melt

ggplot(data_dist, aes(x= value)) + 
  geom_density(fill= "#2b8cbe") + 
  facet_wrap(~variable, scales= "free") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Zmienne, poza length, nie mają rozkładu normalnego.

Korelacja zmiennych

heatmaply(hetcor(as.data.frame(w_duplicates)), k_col = 2, k_row = 3)

Z powodu różnych klas kolumn, np. xmonth jest zmienną kategoryczną, length ciągłą a X porządkową. Została wyliczona heterogeniczna macierz korelacji, metodą hetcor z biblioteki ploycor.

Wykres po czasie

p <- ggplot(w_duplicates, aes(x= X, y= length)) +
  geom_line(alpha= 0.3) +
  geom_smooth(method= "gam", formula= y ~ s(x, k= 100), size= 1) +
  ggtitle("Zmiana rozmiaru złowionego śledzia w czasie")

ggplotly(p)

Jako, że dane zostały uporządkowane chronologicznie, długość śledzia w czasie prezentowany jest przy użyciu liczby porządkowej X. Z powodu ilości danych, który znacznie obniża czytelność wykresu, została zastosowana metoda smooth, która pozwoli na odkrycie ogólnego wzorca. Użycie wygładzenia liniowego, nie byłoby dostatecznie odpowiednie dla zebranego zestawu danych, dlatego został użyty uogólniony model addytywny gam.

Największa korelacja dotyczy par opisujących planktony, lcop1 i chel1 oraz lcop2 i chel2. Duży współczynnik korelacji pomiędzy cumf oraz totaln, przez co możemy wnioskować, iż wraz ze wzrostem łącznej liczby ryb złowionych w połowie rośnie natężenie połowów. Co więcej możemy zaobreswować korelację pomiędzy cumf oraz fbar - łączne roczne natężenie połowów było wysokie tak samo jak ich intensywność.

Regresor

Regresor ma za zadanie przewidzieć rozmiary śledzia w kolejnych połowach. Dane zostały podzielone na dwa zbiory: uczący i testowy, z czego 75% całego zbioru zostało potraktowane jako uczące. Uczenie odbyło się przy użyciu metody Repeated Cross-Validation, z powodu niewielkich różnic wartości w zbiorze zastosowano liczbę powtórzeń na poziomie 5 z liczbą powtórzen 2. Model jest tworzony w opariu o model klasyfikacyjny Random Forrest.

fit <- lm(length ~ ., data = no_x)

# Miara R^2
summary(fit)$r.squared 
## [1] 0.3315561
# Błąd średnio-kwadratowy
rmse <- function(num) sqrt(sum(num^2)/length(num))
rmse(fit$residuals)
## [1] 1.351338
inTraining <- 
    createDataPartition(
        y = no_x$length,
        p = 0.75,
        list = FALSE)

training <- no_x[ inTraining,]
testing  <- no_x[-inTraining,]

ctrl <- trainControl(
    method = "repeatedcv",
    number = 2,
    repeats = 5)

fit <- train(length ~ .,
             data = training,
             method = "rf",
             trControl = ctrl,
             ntree = 2)
## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:ggplot2':
## 
##     margin
fit
## Random Forest 
## 
## 39438 samples
##    14 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 19719, 19719, 19720, 19718, 19718, 19720, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared 
##    2    1.179349  0.4903080
##   13    1.155392  0.5115404
##   24    1.157489  0.5099190
## 
## RMSE was used to select the optimal model using  the smallest value.
## The final value used for the model was mtry = 13.
plot(fit)

rfClasses <- predict(fit, newdata = testing)
summary(rfClasses)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.09   24.56   25.35   25.31   26.24   28.57
df<-data.frame(rfClasses)
ggplot(df, aes_string(x = rfClasses)) + 
  geom_histogram(bins= 100, fill= "#0087BD")  + 
  ggtitle("Przewidywany rozmiar śledzia") + 
  theme_bw() + 
  labs(x= "Rozmiar śledzi", y="Liczba")

Analiza waznosci atrybutów

fit_rf <- randomForest(length ~ ., no_x)
fit_rf
## 
## Call:
##  randomForest(formula = length ~ ., data = no_x) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 4
## 
##           Mean of squared residuals: 1.319012
##                     % Var explained: 51.72
importance_df <- importance(fit_rf)
importance_df <- data.frame(var = rownames(importance_df), importance = importance_df[, 1])
importance_df$var <- factor(importance_df$var, levels = importance_df[order(importance_df$importance), "var"])

ggplot(importance_df, aes(x = var, y = importance)) +
  geom_bar(stat = "identity", fill = "#2b8cbe") + 
  ggtitle("Ważność zmiennych") + 
  theme_bw()

ggplot(data, aes(x = length, y = sst)) + 
  geom_smooth() + 
  ggtitle("Zależność długości śledzia od temperatury przy powierzchni wody") + 
  theme_bw()
## `geom_smooth()` using method = 'gam'

Jak możemy zaobserwować na powyższym wykresie, długość śledzia maleje wraz ze wzrostem temperatury przy powierzchni wody.

ggplot(data, aes(X, sst)) + 
  geom_smooth() + 
  theme_bw() + 
  ggtitle("Zmiana temperatury wody w czasie") + 
  labs(x= "Czas - l.porzadkowa", y="Temperatura[°C]")
## `geom_smooth()` using method = 'gam'

Natomiast temperatura rosła przez ostatnie lata, co spowodowało znaczne obniżenie długości wyławianych śledzi.

ggplot(data, aes(X, nao)) + 
  geom_smooth() + 
  theme_bw() + 
  ggtitle("Zmiana Oscylacji Północnoatlantyckiej w czasie") + 
  labs(x= "Czas - l.porzadkowa", y="Oscylacja Północnoatlantycka")
## `geom_smooth()` using method = 'gam'

Podsumowanie

Konkludując powyższe informacje, możemy postawić diagnozę problemu - w ostatnich latach znacznie wzrosła temperatura przy powierzchni wody, co negatywnie wpłynęło na długość wyławianego śledzia. Wpływ na to ma również zmiana oscylacji północnoatlantyckiej - jest to zjawisko związane z globalną cyrkulacją powietrza i wody oceanicznej, ujawnia się poprzez fluktuacje takich parametrów, jak ciśnienie, temperatura, prędkość wiatru, ilość opadów.